##
# Stacked RDs

### Prep ###

.libPaths( c(.libPaths(), ""))

rm(list=ls())
setwd('')
if (!dir.exists('../2019output/ssn_rd'))   dir.create('../2019output/ssn_rd')
filepath = ''
source("methods.R")


### Helpful methods ###
days_to_base = function(data)  {
   data = data[month(dob) <= 6, base_day := paste0(as.character(year(dob)), "-01-02")]
   data = data[month(dob) >= 7, base_day := paste0(as.character(year(dob) + 1), "-01-02")]
   data = data[, base_day := as.Date(base_day)]
   data = data[, days_since := as.numeric(difftime(dob, base_day, units = 'days'))]
   return(data)
}

normalize_rd = function(rd, ROUND)  {
   rd = rd[days_since == -ROUND, ':='(claim_age_norm = claim_age, retire_age_norm = retire_age)]
   rd = rd[, ':='(claim_age_norm = max(claim_age_norm, na.rm = T), 
                  retire_age_norm = max(retire_age_norm, na.rm = T)), by = 'period']
   rd = rd[, ':='(claim_age_norm = claim_age - claim_age_norm, 
                  retire_age_norm = retire_age - retire_age_norm)]
}
   

### Load data ###

datafile = paste0('../2019data/by_ssn.csv')
dt = data.table( fread(datafile) )
dt = dt[claim_age >= 61.5 & claim_age <= 70]
#dt = dt[claim_age < 61.5 | claim_age > 70, claim_age := NA]     

dt = dt[, dob := as.Date(dob, format = '%d%b%Y')]
dt = days_to_base(dt)
dt = dt[abs(days_since) <= 90 & abs(days_since) >= 2,]

dt = dt[year(base_day) >= 1932 & year(base_day) <= 1943,]
dt = dt[year(base_day) <= 1937, period := '1932-37'][year(base_day) >= 1938, period := '1938-43']
dt = dt[, ':='(	 claimage62 = as.numeric(claim_age >= 62 & claim_age < 63), 
			 claimage63 = as.numeric(claim_age >= 63 & claim_age < 64), 
			 claimage64 = as.numeric(claim_age >= 64 & claim_age < 65), 
			 claimage65 = as.numeric(claim_age >= 65 & claim_age < 66), 
			 claimage66 = as.numeric(claim_age >= 66 & claim_age < 67), 
			 claimage67 = as.numeric(claim_age >= 67 & claim_age < 68),
			 claimagenew2mo = as.numeric(
				(claim_age >= 65 		& claim_age < 65.16 	& year(base_day) <= 1937) |
				(claim_age >= 65.16 	& claim_age < 65.33 	& year(base_day) == 1938) |
				(claim_age >= 65.33 	& claim_age < 65.5 	& year(base_day) == 1939) |
				(claim_age >= 65.5	& claim_age < 65.66	& year(base_day) == 1940) |
				(claim_age >= 65.66	& claim_age < 65.83 	& year(base_day) == 1941) | 
				(claim_age >= 65.83	& claim_age < 66		& year(base_day) == 1942) |
				(claim_age >= 66		& claim_age < 66.16	& year(base_day) >= 1943) ),
			claimagenewgt = as.numeric(
				(claim_age >= 65 		& year(base_day) <= 1937) |
				(claim_age >= 65.16 	& year(base_day) == 1938) |
				(claim_age >= 65.33 	& year(base_day) == 1939) |
				(claim_age >= 65.5	& year(base_day) == 1940) |
				(claim_age >= 65.66	& year(base_day) == 1941) | 
				(claim_age >= 65.83	& year(base_day) == 1942) |
				(claim_age >= 66		& year(base_day) >= 1943) ),
			 male = as.numeric(sex == "M"))]



# Alternative definition of earnings: Use same calendar year for each side of cutoff #
dt = dt[, ':='(	alt_sern62 = ifelse(days_since > 0, sern62, sern63),
			alt_sern63 = ifelse(days_since > 0, sern63, sern64),
			alt_sern64 = ifelse(days_since > 0, sern64, sern65),
			alt_sern65 = ifelse(days_since > 0, sern65, sern66),
			alt_sern66 = ifelse(days_since > 0, sern66, sern67))]


summary(dt)


### Run RD Regressions ###

rdwc = copy(dt)
rdwc = rdwc[, ':='(treated = as.numeric(year(base_day) >= 1938), post = as.numeric(days_since > 0))]
rdwc = rdwc[, ':='(treated_post = treated*post, base_year = year(base_day))]
rdwc = rdwc[, ':='(triangle = (90 - abs(days_since))/90 )]

rdnc = copy(rdwc)
rdnc = rdnc[treated == 1,]

rdcc = copy(rdwc)
rdcc = rdcc[treated == 0,]


for (var in c("claim_age", "claimage62", "claimage63", "claimage64", "claimage65", "claimage66", "claimagenew2mo", "claimagenewgt", "life_income", "alt_sern62", "alt_sern63", "alt_sern64", "alt_sern65", "alt_sern66")) {

	eval(parse(text = paste0("regnc_" ,var, " = felm( ",var," ~ post + post*days_since | 0 | 0 | days_since, data = rdnc)")))
	eval(parse(text = paste0("regnc_" ,var, "_sexM = felm( ",var," ~ post + post*days_since | 0 | 0 | days_since, data = rdnc[sex == 'M'])")))
	eval(parse(text = paste0("regnc_" ,var, "_sexF = felm( ",var," ~ post + post*days_since | 0 | 0 | days_since, data = rdnc[sex == 'F'])")))
	eval(parse(text = paste0("regnc_" ,var, "_white0 = felm( ",var," ~ post + post*days_since | 0 | 0 | days_since, data = rdnc[white == 0])")))
	eval(parse(text = paste0("regnc_" ,var, "_white1 = felm( ",var," ~ post + post*days_since | 0 | 0 | days_since, data = rdnc[white == 1])")))
	eval(parse(text = paste0("regnc_" ,var, "_high0 = felm( ",var," ~ post + post*days_since | 0 | 0 | days_since, data = rdnc[high_income == 0])")))
	eval(parse(text = paste0("regnc_" ,var, "_high1 = felm( ",var," ~ post + post*days_since | 0 | 0 | days_since, data = rdnc[high_income == 1])")))

	eval(parse(text = paste0("regnc_cov_" ,var, " = felm( ",var," ~ post + post*days_since + sex + white + high_income | 0 | 0 | days_since, data = rdnc)")))
	eval(parse(text = paste0("regnc_cov_" ,var, "_sexM = felm( ",var," ~ post + post*days_since + white + high_income | 0 | 0 | days_since, data = rdnc[sex == 'M'])")))
	eval(parse(text = paste0("regnc_cov_" ,var, "_sexF = felm( ",var," ~ post + post*days_since + white + high_income | 0 | 0 | days_since, data = rdnc[sex == 'F'])")))
	eval(parse(text = paste0("regnc_cov_" ,var, "_white0 = felm( ",var," ~ post + post*days_since + sex + high_income | 0 | 0 | days_since, data = rdnc[white == 0])")))
	eval(parse(text = paste0("regnc_cov_" ,var, "_white1 = felm( ",var," ~ post + post*days_since + sex + high_income | 0 | 0 | days_since, data = rdnc[white == 1])")))
	eval(parse(text = paste0("regnc_cov_" ,var, "_high0 = felm( ",var," ~ post + post*days_since + sex + white | 0 | 0 | days_since, data = rdnc[high_income == 0])")))
	eval(parse(text = paste0("regnc_cov_" ,var, "_high1 = felm( ",var," ~ post + post*days_since + sex + white | 0 | 0 | days_since, data = rdnc[high_income == 1])")))
	
	eval(parse(text = paste0("stargazer(regnc_",var,", regnc_",var,"_sexM, regnc_",var,"_sexF, regnc_",var,"_white0, regnc_",var,"_white1, regnc_",var,"_high0, regnc_",var,"_high1,
          	title = 'RD on ",var," (treatment cohort only)',  keep.stat = c('n','rsq'),
          	add.lines = list(c('Sample', 'All' ,'Male', 'Female', 'Nonwhite', 'White', 'Low Income', 'High Income')),
          	type = 'text', 
		out = paste0(filepath, '/Table2_",var,".txt'))")))

	eval(parse(text = paste0("stargazer(regnc_cov_",var,", regnc_cov_",var,"_sexM, regnc_cov_",var,"_sexF, regnc_cov_",var,"_white0, regnc_cov_",var,"_white1, regnc_cov_",var,"_high0, regnc_cov_",var,"_high1,
          	title = 'RD on ",var," (treatment cohort only)',  keep.stat = c('n','rsq'),
          	add.lines = list(c('Sample', 'All' ,'Male', 'Female', 'Nonwhite', 'White', 'Low Income', 'High Income')),
          	type = 'text', 
		out = paste0(filepath, '/Table2_cov_",var,".txt'))")))

}




